package fem2;

import math2.FunctionBasisOnRnToR;
import math2.FunctionRnToR;
import math2.LinAlgUtilities;
import math2.PolynomialBasisOnRnToR;
import fem2.enu.MeshPartType;

/**
 * Edge load along two node of a mesh part. It takes set of basis function
 * associated with node from the mesh part to integrate the external force.
 * Therefore, it can be used with Xfem
 * 
 * @author hbui
 * 
 */
public class EdgeLoad extends Load {

	private FunctionRnToR[] lcs;
	int[] ids;

	/**
	 * 
	 * @param mp
	 * @param lcs
	 */
	public EdgeLoad(MeshPart mp, Node n1, Node n2, FunctionRnToR... lcs) {
		super(mp);
		assert lcs.length == mp.countFields() : "incompatible load components and active fields on mesh part";
		this.lcs = lcs;

		int i1 = mp.contains(n1);
		int i2 = mp.contains(n2);
		if (i1 != -1 && i2 != -1) {
			if (i1 != i2) {
				ids = new int[] { i1, i2 };
			} else {
				throw new Error("two nodes are the same");
			}
		} else {
			throw new Error("some of the nodes do not belong to the mesh part");
		}
	}

	@Override
	public double[] computeRe(double t) {
		FunctionBasisOnRnToR basis = getMeshPart().getBasis();
		int n = basis.getDimension();
		int d = getMeshPart().countFields();
		Node[] nodes = getMeshPart().getNodes();
		int dim = getMeshPart().getDimension();
		double[] re = new double[n * d];

		/*
		 * select the maximum order for integration
		 */
		int p;
		if (basis instanceof PolynomialBasisOnRnToR) {
			p = ((PolynomialBasisOnRnToR) basis).getP();
		} else {
			p = 10; // set a large p to account for accuracy
		}
		int pn = 0;
		for (int j = 0; j < d; j++) {
			if (this.lcs[j] != null) {
				if (lcs[j].getP() > pn) {
					pn = lcs[j].getP();
				}
			}
		}
		p += pn;

		/*
		 * integrate for basis function associated with each node
		 */
		IntegrationManager integrator = IntegrationManager.getInstance();
		MeshPart mp = new Edge(nodes[ids[0]], nodes[ids[1]]);
		mp.setNumberOfFields(dim);
		MeshPartType type = mp.getType();

		double[] xi1 = IntegrationManager.getInstance().getReferenceCoordinate(
				getMeshPart().getType(), ids[0]);
		double[] xi2 = IntegrationManager.getInstance().getReferenceCoordinate(
				getMeshPart().getType(), ids[1]);

		for (int i = 0; i < ids.length; i++) {
			int cnt = 0;
			for (int j = 0; j < ids[i]; j++) {
				cnt += nodes[j].countModes();
			}

			int nmodes = nodes[ids[i]].countModes();

			for (int j = 0; j < nmodes; j++) {
				int start = cnt + j;
				FunctionRnToR Nj = basis.getBasisFunction(start);

				for (int k = 0; k < d; k++) {
					if (lcs[k] != null) {
						int ng = integrator.getNumIntegrationPoints(type, p);
						for (int l = 0; l < ng; l++) {
							double[] xi = integrator.getIntegrationPoint(type, p, l);
							double w = integrator.getIntegrationWeight(type, p, l);
							double[][] jac = mp.jacXAt(xi);
							double detJac = Math.sqrt(LinAlgUtilities.determinantJJt(jac));

							double[] ne = mp.valueAt(xi);

							double[] xiRef = new double[] { ne[0] * xi1[0] + ne[1] * xi2[0],
									ne[0] * xi1[1] + ne[1] * xi2[1] };

							re[start * d + k] += lcs[k].valueAt(xi) * Nj.valueAt(xiRef)
									* (w * detJac) * lambda.valueAt(t);
						}
					} else {
						re[start * d + k] = 0;
					}
				}
			}

		}

		Debugger.watch(re);
		Debugger.watch("", getMeshPart().getDofNumbers());

		return re;
	}
}
